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ABSTRACT 



Context. Numerical simulations predict that a considerable fraction of the missing baryons at redshift z ~ rest in the so-called warm- 
hot intergalactic medium (WHIM). The filaments and sheets of the WHIM have high temperatures (10 5 - 10 7 K) and a high degree 
of ionization but only low to intermediate densities. Therefore, their reliable detection is a challenging task for today's observational 
cosmology. The particular physical conditions of the WHIM structures, e.g. density and temperature profiles, or velocity fields, are 
expected to leave their special imprint on spectroscopic observations. 

Aims. In order to get further insight into these conditions, we performed hydrodynamical simulations of the WHIM. Instead of 
analyzing extensive simulations of cosmological structure formation, we simulate certain well-defined structures and studied the 
impact of different physical processes as well as of the scale dependencies. 

Methods. We started with a comprehensive study of the one-dimensional collapse (pancake) and examined the influence of radiative 
cooling, heating due to an UV background, and thermal conduction. We investigated the effect of small-scale perturbations given by 
the cosmological power spectrum. 

Results. If the initial perturbation length scale L exceeds « 2 Mpc the collapse leads to shock-confined structures. As a result of 
radiative cooling and of heating due to an UV background a relatively cold and dense core forms in the one-dimensional case. The 
properties of the core (extension, density, and temperature) are correlated with L. For longer L the core sizes are more concentrated. 
Thermal conduction enhances this trend and may even result in an evaporation of the core. Our estimates predict that a core may start 
to evaporate for perturbation lengths longer than L ss 30 Mpc. Though the physics in the corresponding three-dimensional case is 
much more complex, one might expect a similar regulation mechanism with respect to the cold streams along filaments, too. However, 
this question will be addressed in a forthcoming paper. The obtained detailed profiles for density and temperature for prototype WHIM 
structures allow for the determination of possible spectral signatures by the WHIM. 

Key words, cosmology: theory - methods: numerical - hydrodynamics - intergalactic medium 



1. Introduction 

At high redshifts z > 2 most of the baryons in the Universe rest in 
the intergalactic medium (IGM) and can be uniquely described 
as gas at still low-density contrast 5 < 10. It is almost iden- 
tically distributed as the underlying dark matter and is highly 
ionized by the UV background radiation (T < 10 5 K). The sub- 
sequent evolution changes that picture. At redshift z ~ 0, only a 
fraction of as 30% of the IGM i s still existing unde r conditions 
comparable with those at z > 2 dStocke et al. 12004b . During the 
evolution toward low redshifts, the mean-scale streaming mo- 
tions could have led to shock-confined filaments containing gas 
at much higher temperatures. 

Numerical simulations by ICen & Ostrikerl (119991) suggest 
that approximately 30 % to 50 % of the cosmic baryons at 
z — are in the form of the intergalactic medium with a tem- 
perature of 10 5 K < T < 10 7 K, which is called warm-hot inter- 
galactic medium ( WHIM). Further n umerical simulations (e.g., 
iDave et alj 1200 U iDolag et alj 120061) with different numerical 
schemes and resolutions also consistently support this picture. 
These numerical predictions have initiated much observational 
effort in order to reveal the existence of the WHIM. 

Owing to the high degree of ionization, the observational 
signature of the WHIM is very weak, in particular with re- 
spect to neutral hydrogen. Therefore, the detection of highly 



ionized metal lines is much more promising. Observationally. 
the WHIM was first proposed through its metal absorption fea 



hires i n the spectra of bright quasars and blaze rs (Hellsten et al 

: & Canizares 2000; Cen et al 



2001 



1998; Perna & Loeb 1998 



Fang & Bryan 1200 ll) . After the first detection of O vi ab- 



sorptio n li nes in the spectra of a bright quasar by Tripp et aTl 
(2000) and lTrippeta"iT(l200lh . a number of detections were re- 
ported thr ough absorption features of Ovi, Ovn, Ovm, and 
Ne ix ions dNicastro et alJl2002t iFang et al]l2002t iMathur et all 
12003b iFuiimoto et alJl2004f) 7 but they are considered to be rather 
tentative. A detecti on with suffi c iently hi gh signal-to-noise ra- 
tio is reported bv iNicastro et al.l d2005albl) . They found absorp- 
tion signatures of the WHIM at two redshifts in the spectra of 
the blazar Mrk421 during its two outburst phases. Future pro- 
posed missions such as International X- Ray observatory (IXO) 
are expected to detect numerous WHIM absorbers. Detection of 
WHIM absorption in the spectra of a fterglows of gamm a-ray 
bursts (GRBs) were also proposed by Elvis et al. (2004) using 
dedicated m issions such as Pharos and were considered more 
recently by iBranchini et alj d2009l) for the prospects opened 
by the recently propo sed satellite missions EDGE and XENIA. 
iKawahara et alJd2006l) investigated the feasibility of these detec- 
tions in a realistic manner based on cosmological hydrodynamic 
simulations. 
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Additionally, several tentative detections o f the W HIM 
throug h its metal line emissi o n are claimed by iKaastra et al.l 
d2003h and iFinoguenov et al.l d2003l) with the XMM-Newton 
satellite. However, these detections are not significant enough 
to exclude the possibility that the observed emission lines are 
of Galactic origin because of the limited e nergy resolution 
(^ 80 eV) of the current X- ray detectors. lYoshikawa et al.l 
(120031) . lYoshikawa eTal] d2004 . and iFang et al.l (120051) showed 
that future X-ray missions equipped with a high-energy resolu- 
tion spectrograph such as DIOS (Diffuse Intergalactic Oxygen 
Surveyor) and MBE (Missing Baryon Explorer) can convinc- 
ingly detect the line emission of the W HIM. 

A comprehensive review is given by Prochaska & Tumli nsonl 

(2008) . 

It is still an open question how much the WHIM contributes 
to the anisotropics of the cosmic microwave background radia- 
tion via the Sunyaev-Zeldovich effect (SZ-effect). Although the 
density contrast of the WHIM is moderate (6 < 100), its tem- 
perature is high (10 5 K<r<10 7 K), and it is supposed to 
make a significant contribution to the cosmic baryon budget 
of ^ 5 0%. Estimates provided by lAtrio- Barandela & Miicket 
(2006k lAtrio-Barandela et al l (120081) . and iGenova-Santos et al. 

(2009) indicate on a non-negligible contribution, which under 
certain conditions might be even comparable with the over- 
all SZ contribution of clusters of galaxies. Thus the SZ effect 
could serve as an additional detection channel for the WHIM. 
However, the strength of a thermal SZ effect is still a matter of 
debate because the results obtained by numerical simulations are 
much less pronounced. 

Therefore, it is of principal importance to investigate the de- 
tailed thermodynamic state and the internal kinematics of the 
structures which may hide a large fraction of cosmic baryons. 
The detailed physics of the WHIM is highly demanding for com- 
putational astrophysics, however. The treatment of low-density 
regions in great detail is difficult. Compared with the numer- 
ical handling of high-density matter distributions where adap- 
tive techniques can be applied, for low-density regions, nec- 
essary higher overall particle and/or grid number is unavoid- 
able for an appropriate description. In addition, higher resolu- 
tion calls for a more detailed consideration of local physics, e.g., 
star formation, feedback, contamination by heavy elements, etc. 
Altogether the computational effort is immensely more compli- 
cated if considered within a cosmological context ( Schave et al. 
2009; iBertone et al.ll2009h . Still, an adequate treatment of low- 
density regions is despite of the currently available highly de- 
veloped computational techniques at the limit or still beyond the 
near future potential. 

In this paper, we consider the formation of WHIM structures 
at z = 0. Contrary to the extensive and complex treatment within 
the context of cosmological structure formation, we investigate 
the evolution of a single one-dimensional prototype sheet in- 
cluding most of the relevant physical processes. We will deter- 
mine the detailed temperature and density profiles and determine 
which processes the latter might depend on. In a forthcoming 
paper we will extend this study toward three-dimensional struc- 
tures. We will use a similar approach as presented here to rep- 
resent a halo - filament scenario. This will enable us to obtain 
further information, in particular about the spatial velocity field. 

Although we avoided performing very time-consuming full 
cosmological simulations we obtain sufficient reliable informa- 
tion about the pre- and post-shock evolution of the WHIM struc- 
tures depending on the assumed initial conditions (amplitude and 
scale-size of initial perturbations). The temperature and density 
profiles are mainly defined by the latter parameters and the phys- 



ical processes that are incorporated. The knowledge of the tem- 
perature and density distribution within WHIM sheets is impor- 
tant for the estimate of the probability to observe particular fea- 
tures (tracers) in the related spectra. This concerns the principal 
observability of certain lines or combinations of lines as well as 
the probability for observations (covering factor) in general. 

The paper is organized as follows: In the next section we 
give the rationale for the considerations given in the paper and 
address our basic assumptions. In Sect. [3] we describe the sys- 
tem of our hydrodynamical equations and basic relations. In 
Appendix lAl we provide the system of chemical rate equations 
and the approximated model for the UV background evolution. 
We also shortly describe the numerical code evora, which we 
specifically developed for this study. The details are given in 
Appendix IB1 In Sect. [4] we investigate the one-dimensional col- 
lapse (pancake formation) for different incorporated physics in 
detail. The dependence of the results on the initial perturbation 
scale is shown in Sect. [5] We discuss our results and their possi- 
ble implications in the final section. 



2. Pancake formation as a model for the WHIM 

Most of the IGM gas distribution is the result of one-and two- 
dimensional collapse processes. The basic theory for the forma- 
tion and evolution of structur e is already well understood sinc e 
the sixties of the last century (Doroshk evich & Zeldovichl 19641) . 
According to these theories, the most probable formation pro- 
cess starts first with the one-dimensional collapse. Only if the 
underlying dark matter-distribution enters the first caustic, multi- 
streaming of matter leads to the formation two-dimensional fila- 
ments and eventually knots, which are characterized by matter 
collapsing in three dimensions. This evolution is closely fol- 
lowed by the gas distribution. A basic description of the gas 
physics in cluding shock appearance w as given in the pioneering 
work of ISunvaev & Zeldovichl (1 19721) in the context of galaxy 
formation. The directions (orientations) for the one-dimensional 
collapsed sheets are determined by the highest eigenvalue of the 
deformation tensor, which can be attributed to the initial linear 
densit y perturbations. According to Doroshkevich & Shandarin 
(1 19781) the probability that two or even three of the initial eigen- 
values are identical or nearly equal is extremely low. Therefore, 
the one-dimensional collapse, and at a certain evolutionary stage 
the two-dimensional collapse, is the dominating structural evo- 
lution process. 

The WHIM structures under consideration (sheet or fila- 
mentary structure) are supposed to reach the non-linear stage 
of evolution not later than at z = 0. If the perturbation scale 
is large enough, the initially perturbed spatial region remains 
Jeans-unstable throughout the whole cosmological evolution, 
i.e., when the mean IGM temperature was raised to about T » 
10 4 K during the cosmic reionization. In order to form shocks, 
the infall velocity of the collapsing gas must reach the speed of 
sound or even go beyond. These conditions lead us to pertur- 
bations on scales initially larger than 1-2 Mpc comoving. The 
structures arising from those perturbations are expected to form 
the large-scale network of the WHIM. 

We start our investigation with the consideration of the one- 
dimensional collapse of one perturbation with a given length 
scale. This scenario is commonly referred to as cosmic pan- 
cake formation. The particular importance of st udying the one- 
dimen sional planar collapse was stressed by Struck-M arcelll 
( 1988). The advantages to restrict ourselves to the simplest geo- 
metrical structures are obvious: 
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- The one-dimensional pancake formation describes the most 
common collapse formation process in the universe at large 
scales (> 1 Mpc). 

- It describes the preliminary phase (predecessor) of collapse 
processes of higher dimensions. 

- It allows a spatial resolution far beyond recent capabilities 
for three-dimensional simulations. 

- For special cases analytical so lutions are availab l e. Those 
may serve as test cases (e.g., iBrvan etall Il995t iTevssiej 
120021) . 

- The high symmetry of the considered configurations does 
not influence the physical state and the principal distributions 
of temperature and density, but it allows for an additional 
check for possible numerical instabilities and deviations of 
non-physical origin. 

- It allows us to investigate and to control the influence of var- 
ious, subsequently introduced, energetic processes onto the 
thermodynamical evolution. 

We include all relevant processes of radiation cooling as well 
as the heating by photoionization due to the evolving cosmic UV 
background. To that purpose we self-consistently compute the 
abundances of the different ionization levels assuming a primor- 
dial medium, i.e. consisting of hydrogen (H) and helium (He) 
only. Besides using the common assumption of ionization equi- 
librium (IE), we perform simulations where the chemical com- 
position is directly computed by integration of the correspond- 
ing differential equations, and check for differences between the 
two cases. While in many cases IE is a valid assumption, de- 
viations from equilibrium may occur, in particular at low den- 
sities. Furthermore we consider the rol e of thermal condu ction. 
Heat conduction was also included by iBond et al. (1984) con- 
sidering the problem of cooling pancakes with analytical meth- 
ods. At certain conditions very high temperature gradients are 
expected to occur. This may happen particularly in the vicinity 
of shock fronts. Then, thermal conduction can lead to a consid- 
erable change of the temperature profiles. In the cases to be con- 
sidered here, the gas is almost fully ioni zed and t he exp ression 
for the heat conduction coefficient from ISarazinl (1 19881) can be 
used. 

The distribution of the gas in the structures of the WHIM 
(sheets and filaments) is supposed to be very close to that of the 
dark matter. This is true even at late evolution stages. In order to 
simplify the calculations, we consider the baryonic content of the 
Universe only and therefore decrease the number of equations to 
be solved. Because this would neglect the gravitational mass of 
the dark matter, we assume that the dark matter obeys the same 
spatial distribution as the baryons. Concordantly, we rescale the 
baryonic density to the total cosmic matter density when com- 
puting the gravitational potential. For test cases, we checked our 
results for deviations from the solutions including the full dark 
matter dynami cs. For that purpose, we used the RAMSES code 
(Tevssier 2002) and appropriate initial conditions. The deviation 
is negligible for the structures considered here, 

Though considering preferentially one-mode perturbations, 
we also investigate up to which degree small-scale perturbations 
may affect the results. For that purpose, we add Gaussian ran- 
dom perturbations according to the cosmological initial power- 
spectrum. The spatial scale sizes of these fluctuations are lower 
compared to that of the considered large-scale single mode, but 
much higher than the (comoving) initial Jeans length immedi- 
ately after reionization. We will compare the various resulting 
density and temperature profiles. 



3. Theoretical framework 

3.1. Dynamical equations 

We use the standard approach for describing the baryonic com- 
ponent of the universe. The assumed ideal polytropic fluid is 
described by the Euler equations, which may be considered as 
conservation laws for the conserved quantities: the density p, 
the momentum densities pu (u denotes the vector of velocities) 
and the energy density E. The latter is the sum of the kinetic 
energy density £jy„ = l/2p|u| 2 and the internal energy density 
E t h. The internal energy density is related to the pressure p by 
the polytropic equation of state p = (y — 1) E t h, where y denotes 
the adiabatic index of the gas. Throughout this paper we use the 
adiabatic coefficient for a mono-atomic gas, i.e., y = 5/3. 

In order to include the cosmological expa nsion into our sim- 
ulatio ns we use supercomoving coordinates (iMartel & Sh apiro 
1998). This is a transformation of the physical coordinate r and 
the time t into the supercomoving coordinates x = (x, y, z) = r/a 
and the conformal time dr = dt/a 2 , where a is the cosmologi- 
cal expansion factor computed from the Friedman equation. We 
use a ACDM cosmology with the parameters deriv ed from the 
five-year WMAP observations dKomatsu et al.l2009l) £l\ = 0.73, 
Q m = 0.27, and Hq = 71 Mpc km 1 s _1 . The transformed Euler- 



equations that are used throughout this paper are 

^+V-(pu) = (1) 
or 

^£-1. + V • (pu <g> u) + Vp = -pV0 (2) 

OT 

^ + V • (u(E +p)) = -pu • + (F - A) - V • j (3) 

OT 

^.+V-(Su) = Izi( ( r-A)-V-j) . (4) 

OT OV 1 



These equations already include the change in energy owing to 
the heating function F, the cooling function A, and the heat flux 
j caused by thermal conduction. The gravitational potential <p is 
computed using the supercomoving version of Poisson's equa- 
tion 

A4> = AnGa Aot ~^ , (5) 
P 

where p denotes the uniform background density of the Universe 
and p tot the total matter density (baryons + dark matter). As 
already mentioned, we assume similar spatial distributions for 
the dark matter and the baryons, and therefore compute the to- 
tal matter density by p tot = pi fn assuming a baryon fraction of 
f B =£VQ m = 0.16. 

Equation © describes the evolution of the modified entropy 
density S = p/p 7 ~ l , which is necessary for the dual-energy for- 
malism described in the Appendix IB. 21 It can be derived from 
Eq. ([T]-[3]l. The only difference of Eq. (Q]-[5]) with respect to the 
non-comoving equations is the factor a in Eq. © (an additional 
drag term would occur in Eq. (01 and Eq. ((4]) for y + 5/3). An 
extensive derivatio n of the supercomoving coordinates is given 
in the Appendix of iDoumler & Knebel (l2010h . 

In order to follow the chemical network an equation of con- 
tinuity for the number densities n, is needed: 

Qfi ■ 

— ^ + V • (inu) = H, . (6) 
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where H, denotes the source term due to chemical processes. The 
index i indicates the five different species Hi, Hn, Hei, Hen, 
and He m. The electron number density can be computed using 
charge conservation: 



n e - "HII + «Hen + 2« He 



(7) 



Initially, the number densities can be computed from the density 
by rii - XiP/mi, where Xi denotes the primordial mass fraction of 
Hydrogen xh = 0.76 or Helium = 0.24, respectively, and m, 
is the corresponding atomic mass. The temperature necessary for 
computing the different rates of the chemical network is given by 



(8) 



It is often assumed in cosmological simulations, that the time 
scales of the chemical processes are much shorter than the dy- 
namical times. The system is therefore assumed to be in ioniza- 
tion equilibrium (IE). Then, the left-hand side of Eq. © van- 
ishes and the number densities can be computed locally using 
H; = 0. Without that assumption, the chemical source terms have 
to be implemented self-consistently using the full set of Eq. ©. 
This includes the hydrodynamical advection given by the sec- 
ond term. In both cases the chemical source term, as well as 
cooling and heating, are computed from photoionization, colli- 
sional ionization and recombination, dielectric recombination of 
He n, collisional excitation of H i and He n, and bremsstrahlung 
(three-body pr ocesses are n egle cted). The corre sponding rates 
are taken from IB lack) (Il98lh and iKatz et all (Il996l) . The details 
of the chemical network and the model for the UV background 
are presented in AppendixlAl 

Thermal c onduc tion is implemented as presented in 
lJubelgas et al.l d2004h . The heat flux is computed by j = 
-k(T)VT, where k is the heat cond uction coefficient. We use 
the coefficient given in lSarazinl d 1 9881) deri ved from the cl assical 
thermal conductivity due to electrons from Spitzer ( 1962): 



« = 4.6xl0 13 (-J-f (^f erg, 
\10 8 K7 \ 40 J g 



s~' cm" 1 K 



(9) 



where InA is the Coulomb logarithm. We use InA = 37.8. At 
low densities the mean free path of the electrons A e given by 



= 0023 ( T^k) 2 ( To^f^) ' Mpc 



(10) 



can approach the scale length of the temperature gradient At - 
T/\VT\. Then the heat flux becomes saturate d. To take this effect 
into account, we use an effective coefficient (Sarazin 1988) 



*eff = 



1 + 4.2 A e I A T 



(11) 



This approach neglects the influence of magnetic fields. They 
mainly affect dense objects like clusters of galaxies. We focus 
on regions of low to intermediate density. Therefore this approx- 
imation is sufficient for our study. 

3.2. Numerical implementation 

In order to compute the time evolution of the fluid governed by 
the equations introduced in Sect. 13. II we developed the new sim- 
ulation code evora. The code uses a regular grid of equally 
sized cells and evolves cell-averaged quantities over discrete 
time steps. Furthermore, the set of equations is split into four 



subproblems: hydrodynamic advection, gravitational accelera- 
tion, integration of the chemical network, and thermal conduc- 
tion. These problems are solved successively at every time step, 
and every solver uses the quantities updated by its predecessor 
as its input. Here we give a brief overview of the four solvers: 

1 . The hydrodynamic problem is solved using the well-known 
MUSCL scheme in combination with the MINMOD slope 
li mite r and the approximate HLLC Riemann solver (see lTorol 
119991) . 

2. The gravitational potential is computed from Poisson's equa- 
tion using Fourier transfor ms and Green's function (see 
iHocknev & Eastwoodlll988l) . 

3. Under the assumption of IE, the number densities are com- 
puted by iteration from S,- = 0. In the non-IE case they 
are inte grated self-consistentl y using the modified Patankar 
scheme (Burchard et a 1. 2003). In both cases the cooling and 
heating functions are ap plied to the pres sure using the regu- 
lar Patankar scheme (see lPatankarll 19801) . 

4. Thermal conduction is included into our simulations by us- 
ing a central difference scheme. 

A more detailed description of our code can be found in 
Appendix iBl 



4. One-dimensional collapse 

4.1. Gravitohydrodynamics 

In a first step we consider the formation of a pancake without the 
inclusion of any cooling and heating or of thermal conduction. 
This serves as a basic configuration to be compared with the 
results after successively taking into account relevant physical 
processes. The configuration is initialized at z = 99 by imposing 
a single sinusoidal perturbation with an initial amplitude A and 
a wavenumber k = 2n/L (where L denotes the length scale of 
the perturbation) with respect to the background density of the 
Universe. This density distribution is then scaled to the baryonic 
density using fg: 



(12) 



p= -L (A cos(fcx) +p) . 

JB 

According to the linear perturbation theory (see |Peebleslll993[) 
we obtain the corresponding velocity field 



/ d A sin(£ x) 



(13) 



where / = (dlogD)/(dloga), with D denoting the growth fac- 
tor. The initial temperature is set to 100 K. We always choose 
the size of our computational domain to be equal to the pertur- 
bation scale L. For any of our simulations, periodic boundary 
conditions are imposed. 

Before z ~ 1 the medium undergoes adiabatic contraction, 
resulting in a sharp density peak in the center of the box. When 
the local speed of sound matches the infall velocity, two shocks 
form and confine a region of high temperature. In Fig. Q] we show 
density, velocity, pressure, and temperature profiles at different 
redshifts after the moment of shock formation. While moving 
outward, infalling cold gas passes these shocks and its kinetic 
energy is transformed into heat. The associated strong decline in 
velocity is visible in the corresponding profile. As a result, the 
temperature inside the shocked region is several orders of magni- 
tude higher then outside. This process is commonly called shock 
heating. While passing the shock, the gas looses most of its ve- 
locity and does not move further inward, but is accumulated at 
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Fig. 1. Pancake formation without cooling, heating and thermal 
conduction: Profiles for different redshifts. Panel A: Density; 
Panel B: Velocity; Panel C: Pressure; Panel D: Temperature. 
This simulation uses an initial amplitude of A, = 0.02 at z — 99, 
a perturbation scale of L = 8 Mpc, and 16000 grid points. 



the outer wings of the profile, leaving the inner part unchanged. 
With time, a continuously decreasing fraction of matter remains 
outside of the shocked region. The accretion onto the pancake 
declines over time. This results in a slower shock speed and in 
a declining density profile. The final density profile bound by 
the shocks covers about 2.5 orders of magnitud e, and is pro- 
portional to r~ 2//3 (Shan darin & Zeldovichl ll989). The pressure 
profile remains almost constant. This is an expected behavior 
since pressure gradients would be quickly erased by hydrody- 
namic advection. The small deviation from uniformity as well as 
the weak redshift dependence are the imprint of the gravitational 
potential and the cosmological expansion. This almost isobaric 
behavior can be used to explain the shape of the temperature 
profiles. Given a constant pressure, Eq. (0 implies an inverted 
behavior between temperature and density. The temperature is 
given in physical units implying a cosmic evolution oc gt 2 . 

Without heating and cooling, the physical dimensions can be 
eliminated from the hydrodynamic equations, and therefore the 
qualitative outcome of these simulations does not depend on the 
given length scale L. If we impose a reference system given by L, 
the background density of the Universe p, and the Hubble time 



cale - Hi L 2 oc L 2 



(14) 



Thus the temperature in the shocked region scales as the square 
of the length scale of the initial perturbation. 

Besides the length scale of the perturbation, the initial ampli- 
tude is set as a paramet er. Its value d etermines the time of caustic 
formation, as shown in lBrvan et al.l d 1995b . The chosen value of 
A = 0.02 corresponds to a shock formation at redshift z ~ 1, 
which could be a reasonable value for the WHIM. Owing to the 
cosmological expansion, the temperature declines very fast from 
its initial value. Therefore, the initial temperature has a negligi- 
ble influence on the dynamical evolution and on the resulting 
profiles. 

4.2. Cooling and heating 

If cooling and heating are included into the consideration, an in- 
trinsic physical scale is introduced. Unlike before, the physical 
dimension cannot be eliminated from the dynamical equations. 
Simulations using different perturbation scales L differ not only 
quantitatively but also qualitatively. As a consequence, the con- 
straint on the ratio between the s patial resolution Ax an d the local 
Jeans length Aj as presented in lTruelove et alJ (1 19971) has to be 
fulfilled: 



0.25>^=A*J^ 
Aj V nyp 



(15) 



If this criterion is not fulfilled, the pressure is too weak to even 
out small scale perturbations, which are caused by the finite nu- 
merical resolution. These perturbations may then grow rapidly, 
and induce fragmentation for purely numerical reasons. In one- 
dimensional simulations this violates the spatial symmetry of the 
configuration. Without the inclusion of either a (formal) heating 
source or an artificial pressure floor, catastrophic cooling in the 
center will appear, i.e., the density increases while the pressure 
decreases. This will always result in a Jeans length that violates 
Eq. (fT~5T >. For the one-dimensional collapse, the heating due to 
the UV background is sufficient to prevent such a cooling catas- 
trophe. 

In Fig. |2]we present the outcome of our simulations includ- 
ing radiative cooling and heating for different length scales L of 
the initial perturbation. For the computation of the chemical net- 
work, we assume IE, as described in Sect. 13.11 The influence of 
non-IE will be discussed in Sect. 14.41 Like before, the initial am- 
plitude is A — 0.02. With increasing L the number of grid points 
increases from 2000 to 64000, thus keeping a constant spatial 
resolution of 0.5 kpc. Displayed are the density, the density flux 
(instead of the velocity, because it emphasizes the high-density 
region in the center), and the related temperature profiles. We 
choose a logarithmic x-axis, thus focusing on the center of the 
simulation box. 

As an immediate effect of the heating due to the photoion- 
izing UV background, the configuration is heated up to tem- 
peratures of T « 2 x 10 4 K during the reionization at redshift 
z ~ 6. This results in a pressure several orders of magnitudes 
higher than in the non-radiative simulations. Therefore, the adi- 
abatic collapse before redshift (z ~ 1) does not produce one sin- 
gle peak, but an isothermal core supported by pressure. The fur- 
ther evolution now depends on the size of the perturbation length 
scale L. 
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Fig. 2. Pancake formation including cooling and heating for an initial perturbation amplitude of A — 0.02 and for a series of 
perturbation scales of L = (2, 8, 12, 16, 32) Mpc/h comoving. Shown is the outcome of the simulations at z = including cooling 
and heating (solid lines), and without dissipation (gray lines). The profiles are displayed using logarithmic coordinate axes. First 
row: Density profiles; Second row: Density flux profiles; Third row: Temperature profiles. 



For the smallest perturbation scale L — 2 Mpc the speed of 
sound inside this core remains always higher than the infall ve- 
locity, and therefore the shock cannot form anymore. The whole 
profile, now sustained by the pressure of the gas only, is more ex- 
tended than in the non-radiative case. For L > 2 Mpc the infall 
velocity becomes higher than the sound velocity at some mo- 
ment (this can be obtained from the scaling considerations dis- 
cussed above). Thus, a shock is able to form. This shock is not 
generated in the center, but forms at the edges of the pre-shock 
core. From there, it moves outward, like in the non-radiative 
case. Additionally, a fan-like wave penetrates into the core, ef- 
fectively shrinking its size. The whole configuration can be sep- 
arated into an inner isothermal core, a shocked region of higher 
temperatures, and an outer region at low density and low temper- 
ature. The size of the core is decreasing with increasing length 
scale L. Outside of the core region, the results of the simulations 
differ only slightly from the non-radiative case. Especially the 
position of the outer shock appears to be unaffected. A special 
situation occurs in the L — 4 Mpc simulation, where an effec- 
tive outflow can be noticed, which appears as a positive density 
flux outside the core. This is caused by the lower density inside 
the core compared to the runs with higher L, resulting in longer 
cooling times, and thus a less effective dissipation of the energy 
input by the further infall. 



The influence of the perturbation scale on the properties of 
the core will be further examined in Sect. 

4.3. Thermal conduction 

The inclusion of thermal conduction leaves an impact on the evo- 
lution of the pancake only for perturbation scales of L > 30 Mpc. 
The used thermal conduction coefficient k shows a steep temper- 
ature dependence of k oc T 25 . Moreover the temperature scales 
approximately like T oc L 2 . This implies a relatively sharp tran- 
sition between perturbation scales where thermal conduction is 
effective or not. Furthermore, thermal conduction is only impor- 
tant if steep temperature gradients occur. This is particularly ex- 
pected within the immediate neighborhood of the shock fronts 
or/and for the temperature step at the core edges. Then it might 
occur that the thermal heat conduction exceeds the cooling and 
the formed core is heated up (evaporation). For this case, a rough 
order-of-magnitude estimate gives the condition under which 
thermal conduction may overcome the cooling, i.e. -V ■ j > A. 
For an estimate, we assume an average cooling comparable with 
that by recombination or collisional excitation. Then we get the 
relation 

neAj K 10 ' 9 ( To*) cnr2 ' (16) 
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Fig. 3. Influence of thermal conduction is shown at regions with 
suficciently steep temperature gradients (L = 32 Mpc). The gray 
lines show the results of the corresponding simulations without 
thermal conduction. Panel A andB: The density and temperature 
profiles at the outer shock fronts are shifted with respect to the 
case without heat conduction. The high-temperature region is 
extended outwards. The density is increased but over a lower 
volume (the density shock is shifted inwards). Panel C and D: 
Heat conduction leads to a smaller core. 



where n e , Aj are the electron number density and the character- 
istic scale for the temperature gradient, respectively. Owing to 
the considerable temperature difference throughout the transi- 
tion zone, one should take T at the at the core edge, whereas for 
n e one should take the value inside the core. At temperature dif- 
ferences of about 10 6 K and n e « 10~ 4 cirT 3 we get a few kpc 
for the transition scale Aj. 

In Fig.|3]we show the details of the density and temperature 
profiles for one perturbation scale L — 32 Mpc. The influence 
of L will be further discussed in Sect. [5] The top panels show 
the region of the confining shock. In the simulations including 
thermal conduction, the shocks in density and temperature do 
not coincide anymore. Thermal energy from the shocked region 
has been transported outward, heating the domain in front of the 
shock to temperatures comparable to the shocked region. This 
energy is lost by the shocked region, resulting in a lower pressure 
and causing a more confined density profile. The now higher 
temperature in front of the shock results in a higher pressure 
there, slowing down the infalling gas. This deceleration causes 
the small increase in density adjacent to the shock. 

Thermal conduction also affects the resulting core profile in 
the simulations including heating and cooling (see Sect. I4.2l i. 
The inner part of the pancake profile is shown in the bottom pan- 
els of Fig. |3] Now, the core is distinguished by a steep increase 
of the temperature at the core border of approximately one or- 
der in magnitude. There, thermal conduction transports energy 
toward the center of the profile. This additional energy raises 
the pressure in the center, causing an outflow. This results in a 
smaller core size with respect to the simulations without heat 



Fig. 4. Top panel: Chemical composition of a simulated pancake 
(L = 32 Mpc) including cooling and heating without thermal 
conduction (solid lines) and with thermal conduction (dashed 
lines). The colors indicate different species (lines from top to 
bottom): red: H II, cyan: He III, magenta: He II, green: H I, blue: 
He I; Bottom Panel: Chemical composition in the region of the 
shock (indicated by the box in the top panel) using non-IE chem- 
istry for HI (left) and Hell (right). The corresponding profiles 
from the top panel are shown in gray for comparison. [See the 
electronic edition of the Journal for a color version of this fig- 
ure.] 



conduction. Contrary to the outer shock front, there is no dis- 
placement between the density and the temperature profiles. The 
core shrinks as a whole. 

4.4. Chemical composition 

Under quite general conditions, the characteristic time scales for 
chemical processes, e.g., such as ionization and recombination, 
are much shorter than the dynamical ones. Then, it is entirely 
justified to assume IE. The chemical composition is entirely de- 
termined by the temperature and the number densities of species 
engaged within the processes. In the top panel of Fig.|4]we show 
the number density profiles assuming IE for a perturbation scale 
of L = 32 Mpc (this corresponds to the results shown in the 
fourth column in Fig. |2j. The gas is almost completely ionized, 
which leads to very low abundances of H i, He i, and He n. The 
shapes of the profiles resemble the density profile, except for a 
noticeable step at the position of the shock caused by the steep 
decrease in temperature at this position. In the simulations in- 
cluding thermal conduction, the behavior at the position of the 
shock is even more complex: The number densities of the H i, 
He i, and Hen show a gap in the region adjacent to the shock. 
Because of the offset between the temperature shock and the 
density shock caused by thermal conduction (cp. with the pre- 
ceding section), this region combines a low density with a high 
temperature, which causes a higher degree of ionization. 
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Under certain conditions the characteristic time scales be- 
come comparable and the assumption of ionization equilibrium 
may become inappropriate. This might happen at very small 
number densities. In this case, we have to follow the detailed 
evolution of the chemical network using the full set of non-IE 
equations ©. 

In those simulations, omitting the assumption of IE, the re- 
sults differ only marginally from that ones of the IE-simulations. 
Effects on the hydrodynamic evolution, coupled by the cool- 
ing/heating function to the chemical network, cannot be detected 
at all. However, the chemical composition shows a slight devi- 
ation with respect to the IE simulations, and this occurs only in 
the direct vicinity of the shock. In the bottom panels of Fig. [4] 
we show the number density of H i and He n around the shock 
using IE and non-IE. In the non-IE simulation a very small re- 
gion adjacent to the shock exists where for Hen the degree of 
ionization is lower than in the IE case. The corresponding chem- 
ical timescale becomes comparable with the hydrodynamical 
timescale. In combination with the motion of the shock, this pro- 
duces the observed behavior. Since, at the temperatures present 
in that region, the chemical rates for the H i are higher, the chem- 
ical time-scales are short, and a similar feature is not observed. 
Owing to the discussed offset between temperature shock and 
density shock, the chemical time scales at the shock region are 
even larger, which results in a more extended region of delayed 
ionization. The same behavior can be observed for He i, only at 
much lower number densities. 

Though the effects of non-IE for the primordial composition 
are only marginal, the influence of the particular conditions at 
the shock position has been demonstrated. Thus for a medium 
containing some fraction of heavy elements which have very low 
densities for usual abundances, the non-IE must be taken into 
consideration. 

If omitting the assumption of uniform temperature for all 
fluid compo nents (electrons, ions ) the effects on non-IE maybe 
much larger dTevssier et al.lll998l) . 

4.5. Influence of small-scale perturbations 

The initial conditions for cosmological simulations of large- 
scale structure formation are given by a spectrum of perturba- 
tions. To take this into account, we add Gaussian random per- 
turbations according to the cosmological power spectrum to the 
particular perturbation given by Eq. (1121 1. Then, in a given spatial 
region of size L, a pronounced pancake structure will only form 
if the considered initial perturbation amplitude dominates over 
the neighboring perturbation amplitudes at comparable scale 
size. Thus, we subsequently include all perturbation modes up 
to a scale size of (1/8, 1 /4, 1 /2) x L. 

In Fig. [5] we show the density and temperature profiles of 
simulations including the small scale perturbations and cooling 
and heating at z = 0. Clearly, the perturbation modes at scales 
comparable with L have the most impact on the final profiles at 
z — 0. In this case most of the power from perturbation modes 
at neighboring scales will be added. In particular, this leads to 
an enhancement of the core density. The modes smaller than the 
actual Jeans length are erased by the heating due to the UV back- 
ground. However, the magnitudes of the temperatures and their 
coarse profiles in the post-shock region are on the same order 
and the density profiles are nearly preserved. 

We conclude that it is sufficient to consider only the collapse 
of a single (large enough) mode in order to gain coarse infor- 
mation on the thermal and chemical state of the structures of 
interest. In addition, the conservation of the system's symmetry 




-0.2 -0.1 0.1 0.2 -0.2 -0.1 0.1 0.2 

log(x) [Mpc/h] log(x) [Mpc/h] 



Fig. 5. Pancake formation if Gaussian random perturbations are 
taken into account for the initial density field according to the 
cosmological perturbation power spectrum. The resulting den- 
sity (left) and temperature (right) profiles are shown at z — for 
an initial perturbation at L = 8 Mpc. The included upper scale 
length for the Gaussian random perturbations as fraction of L de- 
creases from top to bottom: First row: 1/8 L; Second row: 1/4 L; 
Third row: 1/2 L . The obtained profiles are compared with the 
single-mode consideration given in Fig. [2] The latter is shown in 
gray. 



serves as a probe for the quality of the numerical treatment. The 
non-radiative simulations reproduce the earlier found a nalytic 
results with high accuracy (iShandarin & Zeldovichll989h . 

5. Scaling relations 

In the preceding section it was shown that the length scale of 
the perturbation L is the determining parameter for the evolution 
of the pancake characteristics (temperature and density profiles). 
In Fig. [6] we present the L-dependence of the final values of the 
central temperature T c , the central hydrogen number density «Ho 
the radius of the isothermal core A , and the temperature at its 
edge T . 

For the non-radiative simulations, T c shows the expected be- 
havior: The temperature scales oc L 2 , and thus T c , as well. The 
density does not depend on L, its profile is uniquely oc r~ 2 ^ 3 (see 
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Fig. 6. One-dimensional collapse: Dependence of various quan- 
tities on the perturbation length scale L for simulations without 
(•), with (x) cooling and heating, and with cooling and heat- 
ing and thermal conduction (+). Panel A: Central temperature 
(the dotted line represents T c oc L~° M ; Panel B: Central hy- 
drogen number density (the dotted line represents «h oc L 238 ); 
Panel C: Core size (the dotted line represents A a oc L l 38 ); Panel 
D: Temperature at the core boundary (the dotted line represents 
T oc Z 3 ). In the upper panels the data points with thermal con- 
duction are identical to those without and are therefore omitted. 
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Fig. 8. Top panel: Absolute value of the cooling/heating function 
|r - A| for a hydrogen number density of zih = 4.6 x 1CT 3 cm -3 
corresponding to an overdensity of p/pt = 5000 at z = 0.7 (solid 
line) and for a density of «h = 9.2 x 10~ 4 crrT 3 correspond- 
ing to an overdensity of p/pt, = 1000 (dashed line). The cooling 
function without heating is given in gray for comparison. Bottom 
panel: Equilibrium temperature vs. density. In the range delim- 
ited by the two densities used in the top panel (emphasized by 
the the dashed lines) the relation can be approximated by a power 
lawT„ oc«-° 16 . 
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Fig. 7. Phase space diagram, i.e., the dependence temperature 
versus density, is shown for the quantities at the center of the 
pancake. Shown are curves for L = (16, 24, 32) Mpc using solid, 
dashed, and dotted lines, respectively. The upper gray line shows 
equilibrium temperature at redshift of z — 0.7 (time of shock for- 
mation) and the lower gray line at z — 0. After a period of linear 
growth (a), the perturbation decouples from the cosmic expan- 
sion and contracts adiabatically (b) until shock formation (at the 
maximum of the temperature). Then, after nearly isochoric (c) 
and isobaric (d) evolution stages the central gas arrives at ther- 
mal equilibrium (e). 



Shandarin & Zeldovich|[T989l) . However, because we increase 
the number of grid points in order to keep the spatial resolu- 
tion fixed, larger simulations resolve the central peak better. In 
the result, we get an apparent dependence of the central density, 
i.e., of the density at the innermost grid point, oc L 2 ^ 3 . However, 
the latter reflects only the degree of resolution. 

If including cooling and heating processes, the above rela- 
tions are no longer valid. The central temperature stays roughly 
constant at about 2 x 10 4 K, weakly decreasing at higher L pro- 
portional to Zr - 38 . The central hydrogen number density shows a 
strong dependence on L approaching the asymptote «Hc k L 23& . 
At the edge of the isothermal core, the temperature strongly in- 
creases outwards, but then flattens again. We identify the core 
radius A a as the distance from the center to the inflexion point 
of the temperature profile. Then the core radius shows a depen- 
dence on the scale length proportional to L l 38 . 

The above mentioned scaling relations for T c , «Hc, and 
A a can be explained using simple thermodynamical arguments. 
Actually, the size of the core and its density is fully determined 
by conditions of hydrostatic equilibrium. In this case Euler's 
equations yield for the central pressure p c ~ nn c (p and Poisson's 
equation yields <p ~ nn c Al. Combining those two estimates we 
obtain 

Pc~£n\ c (17) 

Without heating the actual pressure of the cooling gas is not 
able to withstand the gravity forces, and no distinguishable core 
is forming in these simulations. Including heating, the pressure 
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is raised by several orders of magnitude during reionization at 
z~6. 

In Fig. Qwe show the path with respect to the central point 
of the configuration in a phase space diagram, where the tem- 
perature T c is plotted against the hydrogen number density «h c 
(these are physical, not super-comoving quantities). The curve 
enters the shown domain at the center of the bottom of the plot, 
when reionization heats the gas to « 2 x 10 4 K. From there, the 
density decreases at an almost constant temperature according 
to the linear perturbation growth and cosmological expansion. 
Then the gas dynamics decouples from the Hubble flow and fol- 
lows the adiabatic pressure-density relation p c oc . The shock 
forms at the maximum of p c and T c . During the collapse, the 
energy budget of the central region is dominated by the infalling 
matter and not by cooling and heating. Therefore, the scaling re- 
lation for the pressure derived for the non-radiative simulations 
p oc L 2 also holds here. Using this relation (and y = 5/3) we 
obtain scaling relations for the central density and the central 
temperature at the time of shock formation (denoted by the in- 
dex s): 

n Hs oc P \ ly oc L 2/y = L L2 (18) 

T s oc p s /n m oc L 2 - 2 ' 7 = = L & . (19) 

The formation of the two shocks moving outwards changes the 
situation significantly. The energy supply through infall van- 
ishes, and the subsequent evolution is dominated by cooling 
and heating. Nevertheless the pressure inside the shocked re- 
gion, determined by the potential, stays roughly constant, like 
in the case without cooling and heating. Therefore the relation 
p c oc L 2 and Eq. dTTb remains valid. The efficient radiative cool- 
ing decreases the temperature toward the equilibrium tempera- 
ture T e of the cooling/heating function, which is the tempera- 
ture where the contributions of the cooling processes and the 
UV background heating are canceling each other. The gas with a 
higher temperature cools down toward T e , while colder gas gets 
heated. In case of ionizational equilibrium this temperature is a 
function of the density only. In the top panel of Fig. [8] we show 
the absolute value of the normalized cooling/heating function 
|r - A\/n 2 i for number densities corresponding to overdensities 
of 1000 and 5000 at redshift z = 0.7, which is the approximate 
time of shock formation. The equilibrium temperature is given 
by the zero value of \V - A\/n 2 i (the values in Fig.[8]are absolute 
values). A higher density shifts T e toward lower temperatures. 
Computing \F-A\/n 2 i for several densities and tracking the min- 
imum we obtain a relation between «h and T e . This is shown in 
the bottom panel of Fig. [8] This relation is also shown in Fig. [7] 
(solid gray lines), for both z = 0.7 and z = 0. Evidently, when 
reaching T e , the central state remains there, evolving further only 
due to the cosmological expansion. In the interval of overdensi- 
ties of 1000-5000 (at z = 0.7), the relation can be approximated 
by a power law T e oc « H 016 (indicated in Fig. [8). At lower densi- 
ties the behavior deviates only weakly from that relation. Using 
this approximation and p c oc nn c T c oc n^A 2 , oc L 2 we obtain 

«Hc <* Pc/T c oc p e /n H e 16 "He <* L 2M (20) 

T c oc < 16 ocL-°- 38 (21) 

A a oc L/n He <* £~ L38 • (22) 

Although we have only used the approximated relation for the 
equilibrium temperature, the above consideration leads us to ex- 
pressions that fully agree with the results shown in Fig. [6] This 
means that the system eventually approaches a quasi-equilibrium 
state for sufficiently large scales, at least. 



The core size decreases with increasing scale length L ac- 
cording to A a oc L l ■ 38 . As can be seen in Figs.[3]and|6]the heat 
conduction even amplifies this tendency. Thus it might be ex- 
pected that for some scale length L an evaporation of the core 
will happen. Obviously, that will be the case if the core accord- 
ing to the cooling/heating processes becomes of comparable size 
as the above introduced heat transition scale At, i.e. At > A . 
Using the expression ( fT6l l and the derived scaling relations nor- 
malized to the corresponding numerical values at L = 16 Mpc 
we obtain for the scale length at which the core is significantly 
affected by thermal conduction 

L w 30 Mpc . (23) 

At larger scales even evaporation of the core becomes possi- 
ble. This a l so pri ncipally agrees with the results obtained by 
iBond etal.ldl984l) . 

6. Summary and conclusions 

The aim of the present paper is to investigate the detailed physics 
and thermodynamics of a matter distribution which reflects the 
basic properties of structures of the WHIM. From cosmologi- 
cal simulations we know that the WHIM is distributed in sheets 
and filaments. For the description of a matter distribution at low 
density the resolution is a key issue. Therefore, we initially con- 
centrated on the description of the one-dimensional collapse to- 
ward a gaseous sheet at extremely high resolution. Although 
neglecting any interaction between the perturbations on vari- 
ous scales we started with perturbation parameters consistent 
with the cosmological density field. Above an initial perturba- 
tion scale of as 2 Mpc the collapsing gas shocks and further ther- 
modynamics are mainly determined by the cooling and heating 
processes. In case of the one-dimensional collapse the heating 
by the photoionizing UV background is sufficient to keep the 
whole structure in a quasi-equilibrium state. Then the final den- 
sity and temperature distribution depends on the initial pertur- 
bation length scale L only, i.e. all quantities characterizing the 
quasi-equilibrium state may be roughly described as function of 
that length scale L. If keeping the resolution high enough and 
fixed, and increasing the initial length scale L then this goes be- 
yond the numerical capabilities. The obtained scaling relations 
can be used thought to extrapolate the results towards higher L. 

The one-dimensional density and temperature profiles are 
characterized by the existence of a cold and dense core region 
which is at thermal equilibrium at a temperature of T c ~ 2 x 10 4 
K. This core is built even before shock formation, and its proper- 
ties are given by the interplay between radiative cooling and the 
energetic input from the UV background. 

In this respect the question about a possible shielding of the 
UV flux is essential. However, the estimates of the optical depth 
T p h with respect to photoionization by the UV flux show for all 
possible cases r p h <k 1 . The situation may be different for the 
three-dimensional case when higher densities are reached. 

The determining variable is again the perturbation scale L: 
larger collapsing scales lead to a spatially narrower cold region. 
The size of the core region decreases even more if thermal con- 
ductivity becomes efficient. For large enough scales L the tem- 
perature gradients at the transition from the cold core towards 
the shock heated gas are large. In the result, thermal conductiv- 
ity leads to a partial "evaporation" of the core. Using the derived 
scaling relations with respect to the parameter L we can estimate 
the approximate collapse scale when the core will be entirely 
evaporated. In the result, we get a range of scales L between 2 
and 30 Mpc, for which a cold and shock-confined core can exist. 
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Fig. 9. Rendering of the density field at z — for three- 
dimensional simulation using L = 4 Mpc. Overdensities in red, 
underdensities in blue [See the electronic edition of the Journal 
for a color version of this figure.] 



The cold gas in the core region of the WHIM structures must 
contain a relatively large fraction of neutral hydrogen and should 
be detectable as Lya absorption lines in the light of more dis- 
tant sources. If the core size is estimated to be on the order of 
w 10 kpc and the overdensity is in the range 100 < 5 < 1000, 
then Lyff absorption lines could reach column densities of about 
Nm ~ 10 14 to 10 15 cuT 2 . At once, a contribution in absorp- 
tion or emission by the heavy element contamination in the 
WHIM should be detected. A likely d etection of such a sim ulta- 
neous event was reported recently bv lNicastro et alj d2010l) . us- 
ing combined X-ray and optical observations towards the Seyfert 
1 Galaxy PKS 0558-504. Of course, the particular observation 
depends on the actual geometry and orientation of the filament 
and the line of sight. The probability for direct observations of 
the cold core region is rather low. 

Hig h-resolution hydrodynamical simulations o n galaxy for- 
mation (Ocvir k et al.l 2008; IKeres et alj l2009bllal; lAgertz et alj 
2009; Ceverin o et alJl2009l) indicate another very important as- 
pect. At higher redshifts z — 2 - 3, relatively cold gas is flow- 
ing along the filaments toward the connected galactic halos. The 
cold gas streams may even reach inner galactic regions, which 
are closer to the center than the virial shock region. These cold 
gas streams may supply galaxies with fuel for star formation or 
may influence the the galactic star formation, at least. Several 
studies discuss in which way t he existence o f thes e cold streams 
is related to the galaxy mass. IKeres et all d2005l) find that the 
cold gas accretion is only important for t otal galactic masse s 
< 10 11 a Mq. Related estimates are given in iDekel et all d2009h . 
If the WHIM gas is located mainly within a network of large 
shock-confined filaments at z ~ 0, then even recently a cer- 
tain gas flow could be expected towards the corresponding knots 
of the network. Those knots could be matter concentrations of 
galactic or even of galaxy cluster size. 



If extending our above considerations in three dimensions, 
multi-streaming occurs, leading to the formation of a well- 
defined filament and a halo, cp. Fig. [9] Though the shape of the 
gas distribution is extremely idealized, it is expected to exhibit 
the correct density and temperature profiles depending on the 
initial scale length, again. The large-scale gas velocities are de- 
termined by the gravitational potential. In the three-dimensional 
case, the velocities toward the forming filament and even more 
toward the halo are becoming high enough to generate shocks. 
Instead of a core, a relatively cold stream forms along the three- 
dimensional filament, propagating deep into the halo. The tem- 
perature and density profiles for the three-dimensional case are 
much more complex even for the highly idealized geometry. 
We are going to consider this in a forthcoming paper (Klar & 
Miicket, in preparation). However, the existence of a natural reg- 
ulation mechanism for the size and existence of a cold core for 
the one-dimensional case might indicate similar restrictions for 
the cold stream in three dimensions. 
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Appendix A: chemical rates and a model for the UV 
background 

The coefficients used to compute the evolution of the chemical 
network and the cooling and heating function are shown in Table 
I A. 1 1 The values of the r ates of the collisional processes are taken 
from [Katz et afl (Il996l) . w hile the pho toionization and photo- 
heating rates are taken from lBlackl d!981l) . In terms of these rates 
the chemical source term writes: 



Shi = aHii«Hii«« - (fimn e + ym) «hi 
Hhh = (fimn e +ym) n H i - ^hii «hii 

^Hel = (ttHell + ^HeIl)«HeII "e 

- (^Hel«e - yHel) "He I 

^He II = »He III «He III We (/fee I «e ~ THe i) "He I 

- (ttHe II + ^He II) "He II "e 

- (/fee II n e + jHell) "Hell 

^Hein = (/feell"e +7HeIl) "Hell -«HeIIlHHe III "e 



(A.l) 

(A.2) 
(A.3) 



(A.4) 
(A.5) 



The cooling function A is the sum of contributions from the col- 
lisional processes discussed above as well as collisional exci- 
tation of H I and He II and bremsstrahlung, while the heating 
function F is the sum of the heating rates corresponding to the 
photoionization: 

A = £ H I«Hl«e + £HeI«HeI«f + £Hen«HeII«e 

+ rjH II "H II n e + 77He II "He II «f + ^He III "He III «e 



+ W H ii«Hii "« + ^Hl"Hl" C + lAHeII«HeII«e 
+ #("Hn + "HeII +4n H em) «e 
T = £Hl"Hl + £Hel"HeI + £HeII«Hen • 



(A.6) 
(A.7) 



For the computation of photoionization as well as photoheat- 
ing the flux of the UV background is needed. We use a simpli- 
fied model for its redshift dependence which resembles the cur- 
rent view in the literature. ( Gnedi nl2000l:lHaardt & Madaul2001l: 
iBianchi et a l. 2001): 



./() 



10 21 ergs 'cm 2 x{ 



io- 5 

0.5 x 10°- 35 < 6 - d 

1q0.1(3-z) 
1 

0.1 x 10 z 



ifz>8 
if 8 < z < 6 
if 6 <z < 3 
if 3 <z < 1 
if z < 1 



.(A.8) 



A spectrum inversely proportional to the frequency is assumed. 



Appendix B: The evoracode 

Most of the techniques used in our code are common in cos- 
mological simulation codes. We will therefore concentrate on 
specific features of evora. To preserve readability we will give 
one dimensional descriptions in x direction. The generalization 
to three dimensions is straightforward. We will use x' = x(t + At) 
for a quantity x after the time step At and x = x(t) at its begin- 
ning. 



B. 1 . Time stepping 

The length of the time step At is given by the minimum of 
three different constraints. The first is the Courant-Friedrich- 
Levy condition implied by the the hydrodynamic solver: 



At 



cri 



Ax 
\u\ + a 



(B.l) 



where a is the speed of sound. Moreover, the cosmological ex- 
pansion during one single time step is limited by 



a 

At a = — . 
a 



(B.2) 



The last constraint on the time step is associated with thermal 
conduction. Here, the maximal length of the time step can be 
computed by estimating the fraction between the thermal energy 
and and its change due to the thermal conduction: 



Af„ = 



pT pAx 2 



E ± kT/Ax 2 



(B.3) 



Furthermore, each of these time steps is further restricted by 
a heuristic factor < C < 1. We use C c // = 0.5 for the 
CFL-constaint, C a = 0.01 for the cosmological constraint, and 
C tc = 0.9 for the thermal conduction constraint. The time steps 
are computed in every cell, and the minimal value is used to 
compute the global time step. 

The cooling timescale (which is also the chemical timescale) 
can be much shorter than timescales given above. Therefore, it 
is not used as a constraint on the global time step. Instead we use 
subcycling: The time evolution of the number densities and the 
thermal energy density is computed using several shorter time 
steps. This is done locally i n each cell whil e keeping the other 
quantities constant (see also lKav eT al. 2000). 
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Table A.l. Used rates for chemical evolution, cooling, and heating. Rate s of collis ional processes are taken from lKatz et alJ (l996, 
Tables 1 and 2); Photoionization and photoheating rates are taken from Bl ackl (1198 lL Sect. 2.1). T a denotes 77 10" K. 



Process 


Type 


Species 


Coefficient 


Value 




standard recombination 


chemical 


Hn 




8.40 x 10" 


■11 r -0.5 y-0.2 (1 + yO.7-,-1 






Hen 


Q'Hell 


1.50 x 10" 


10 7— 0.6353 






Hem 


Q'Helll 


3.36 x 10" 


-10 r -0.5 r -0.2 (1 + r 0.7 r l 




cooling 


Hn 




8.70 x 10" 


.27 r 0.5r-0.2 (1+r 0.7)-I 
3 ^ 6 






He i 


'/He I 


1.55 x 10" 


-26 y0.3647 






Hen 


'/Hell 


3.48 x 10" 


-26 7-O.5 7— 0.2 /j + jQ.lyl 


dielectric recombination 


chemical 


Hen 




1.90 x 10" 


3 T~ h5 exp(-470000/T)(l + 0.3 exp(-94000/T)) 




cooling 


Hen 


^He II 


1.24 x 10" 


13 r -1 - 5 exp(-470000/T) (1 + 0.3 exp(-94000/T)) 


collisional ionization 


chemical 


Hi 


fain 


5.85 x 10" 


n 7-0.5 + t-o.5^-1 exp(- 157809. 1/r) 






He i 


/^He II 


2.38 x 10" 


n 7-0.5 + 7-0.5^-1 eX p(_286335.4/r) 






Hen 




5.68 x 10" 


12 7-0.5 ^ + 7-0.5-)-] exp(-631515.0/r) 




cooling 


Hi 


<Thi 


1.27 x 10" 


21 7-0.5 (l + T o.Syi eX p(- 157809. 1/r) 






Hen 


^*He II 


1.27 x 10" 


21 7-0.5 ^ + yO.5-,-1 exp (_ 157809. 1/r) 






He m 


^He III 


4.95 x 10" 


22 y-0.5 (J + T 0.5yl eX p(-63 15 15.0/r) 


photoionization 


chemical 


Hi 


7HII 


2.54 x 10" 


jo 






He i 


THell 


2.49 x 10 8 


jo 






Hen 


THem 


1.60 x 10"' 


jo 




heating 


Hi 


£Hi 


7.75 x 10" 


•12 

THi 






He i 


EHel 


2.19 x 10" 


"rne. 






Hen 


£ Hen 


3.10 x 10" 


-11 „ 

7Hell 


collisional excitation 


cooling 


Hi 


<Ahi 


7.50 x 10" 


19 (1 + r 5 - 5 )-' exp(-118348.0/r) 






Hen 




5.54 x 10" 


-17 7-0.397 (1 + r o.5)-i eX p(-473638.0/r) 


bremsstrahlung 


cooling 


all ions 


e 


1.42x10" 


27 8ff T 03 with the gaunt-factor gff = 1.5 



B.2. High Mach-number problem 

In cosmological simulations, one often encounters flows of high 
velocity and low pressure. In these situations, the numerical 
computation of the difference E - E^, needed for the com- 
putation of the pressure, might might not yield reasonable re- 
sults. This is known as the high Mach-number problem. To over- 
come this problem we implement the algorithm suggested by 
Feng et all d2004l) . In addition the conservative quantities we 
also follow the evolution of a modified entropy density S = 
p/p y ~ l . In high Mach flows, where E * E^ n , we use S to com- 
pute the pressure, while elsewhere E is used. After the pressure 
is computed the quantity not used for the computation is recom- 
puted using p, thus keeping both quantities synchronized. 



B.3. Hydrodynamic solver 

If we set the right-hand side of Eq. ([T] - |4j» to zero we obtain the 
homogeneous Euler-equations. These equations are used to com- 
pute the pure hydrodynamic evolution of the fluid. This prob- 
le m is solved u sing the MUSCL-Hancock scheme as presented 
in iTorol {T999). To ensure the mono tonicity of the solution the 
MINMOD slope limiter is applied. The inter-cell fluxes are com- 
puted using a HLLC Riemann-solver. This solver approximates 
the analytic solution of the Riemann problem by three waves 
separating (with velocities S l, S+, Sr) four different states: the 
left (L) and the right (R) initial state, and regions left (*L) and 
right (* R) of the co ntact discontinuity. We use the algorithm 
given in ITorol ( 1 19991 chap. 10.4 and 10.5) and add a prescrip- 
tion for the computation of the modified entropy density and, for 
non-IE simulations, the number densities in the central regions 
★L and *R: 



1i,*K = PK 



> *K 



Pk 



?> K -u K \ p K 

Sk-sJ P y 



S k - uk \ rii,*K 



S k - S * / pk 



(B.5) 



where K = L or K = R for the left or right states. The wave 
speeds are computed by 



S L = min[u L -a L ,u R -a R ] 
Sr = max [ui + aj_, ur + or] 



(B.6) 
(B.7) 



(B.4) 



„ PR- PL+ PLU L (S L~ Ul) ~ PRURiS R~ Ur) 

Pl(S L - U L ) - p R (Sr - Ur) 

Once the state on the interface is determined, the fluxes can be 
computed and the conservative quantities can be updated. 



B.4. Gravitation 

The gravitational potential is obtained though the classical 
algorithm used by Particle-Mesh codes. We solve Poisson's 
equation with Fourier-trans formations and Green's function 
dHocknev & Eastwoodl ll988). Then, the conservative quantities 
are updated with the gravitational source terms. To perform 
the needed Fo urier-transformations w e use the public available 
FFTW library dFrigo & Johnsonll2005l) . 

B.5. Chemical evolution 

The evolution of the number densities and the heating and cool- 
ing originate in the same physical processes and have similar 
timescales. It is therefore necessary to compute their evolution 
in a similar way. In the IE case a solution to H, = can be 
found by iteration. The situation is more difficult in the non-IE 
case. Here, the integration of the system of ordinary differential 
equations h = H, is performed. These are stiff ordinary differ- 
ential equations, and therefore most codes use implicit methods 
for their solution. In our code, we use a different approach and 
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adopt the modified Patankar scheme developed in biochemical to 
oceanography (Burc hard et al.l 120031) . Although it is explicit it . 

ensures the positivity of temperature and number densities and E' = Ei H (j; +1/2 - 7; -1/2) (B.21) 

conserves total amount of hydrogen and helium. The modified ^ 

Euler-Patankar scheme reads: ~ s ; = s . + A* 7_zl (/ . +[/2 _ _ 1/2 ) . (B . 22 ) 



: Hi + Af 



) ,Pik )j d ik — 



At p 



(B.9) 



where p% is the production matrix containing the rates produc- 
ing species i from species k, while dn. is the destruction matrix 
containing the rates that transform species i into species k. That 
implies pn - du and all diagonal coefficients are zero. We ob- 
tain: 

PY2. = aHlI«Hn«e = d 2 \ (B.10) 
P2\ = (jSmne +rHl)«Hl = <^12 (B.ll) 
P34 = (ffHen +^HeIl)nHeII«e = ^43 (B.12) 
P43 = OSHelWc + rHel)"HeI = ^34 (B.13) 
P45 = aHellinHemnf = ^54 (B.14) 
P54 = OSHenWe +yHen)«Hen = ^45- (B.15) 



The other components vanish. If we define 

At p ik x^Atd ik 

P ik = and Di=) . (B.16) 

n k ^ rij 

Eq. dB.9t can be further performed to 

n'. = m + y P ik Xk - Din', = - + Y P ' k x k . (B.17) 

k k 

This is an easy solvable system of linear equations. Since the 
2x2 block matrix for hydrogen and 3x3 matrix for helium are 
not coupled, this system can be solved for them independently. 

In the IE as well as in the non-IE case the up date of the pres - 
sure is performed using original Patankar- Trick (Patankar 1980): 

/ p'\ p + AtT 
p'=p + Af(r-AL) = ^__. (B .18) 



B.6. Thermal conduction 

The algorithm for the thermal conduction is carried out similar to 
the hydrodynamic scheme. After the thermal fluxes between the 
cells are computed, those fluxes are used to update the energy 
density and modified entropy density. Like the hydrodynamic 
scheme this is done in an unsplit fashion. In order to compute 
the thermal flux, first the temperature gradient is computed: 

Ti +i - Ti 

(v,rw= h , (B.i9) 

where i x + 1/2 denotes the position of the interface and i x and 
i x + 1 the cell left and right of the interface. Then, the conduction 
coefficient k and the fraction of mean free path and temperature 
A e /T are computed in cell i x and i x + 1 and then extrapolated to 
the interface i x + 1 /2 by simple averaging. The heat flux is then 

j ix+ i/2 = KV x T (l + ^V.r) ' , (B.20) 

where all quantities are located at i x + 1/2. In a final step, we 
update energy density and modified entropy density according 
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